Naval  Research  Laboratory 

Washington,  DC  20375-5320 


NRL/MR/5341— 16-9680 


A  Parameterized  Pattern-Error 
Objective  for  Large-Scale  Phase  Only 
Array  Pattern  Design 


Dan  P.  Scholnik 

Surveillance  Technology >  Branch 
Radar  Division 


March  21,  2016 


Approved  for  public  release;  distribution  is  unlimited. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  this  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including 
suggestions  for  reducing  this  burden  to  Department  of  Defense,  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports  (0704-0188),  1215  Jefferson  Davis  Highway, 
Suite  1204,  Arlington,  VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  any  penalty  for  failing  to  comply  with  a  collection  of 
information  if  it  does  not  display  a  currently  valid  OMB  control  number.  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1.  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE  3.  DATES  COVERED  (From  -  To) 

2 1  -03-20 1 6  Memorandum  Report  1  October  2012-31  March  20 1 5 


4.  TITLE  AND  SUBTITLE 


5a.  CONTRACT  NUMBER 


A  Parameterized  Pattern-Error  Objective  for  Large-Scale 
Phase-Only  Array  Pattern  Design 


5b.  GRANT  NUMBER 


5c.  PROGRAM  ELEMENT  NUMBER 


6.  AUTHOR(S) 


5d.  PROJECT  NUMBER 


Dan  P.  Scholnik 


5e.  TASK  NUMBER 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 


5f.  WORK  UNIT  NUMBER 

4895 

8.  PERFORMING  ORGANIZATION  REPORT 
NUMBER 


Naval  Research  Laboratory 
4555  Overlook  Avenue,  SW 
Washington,  DC  20375-5320 


NRL/MR/5341— 16-9680 


9.  SPONSORING  /  MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


10.  SPONSOR  /  MONITOR’S  ACRONYM(S) 


Naval  Research  Laboratory 
4555  Overlook  Avenue,  SW 
Washington,  DC  20375-5320 


11.  SPONSOR /MONITOR’S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 


Approved  for  public  release;  distribution  is  unlimited. 


13.  SUPPLEMENTARY  NOTES 


14.  ABSTRACT 

Modern  phased-array  radar  systems  increasingly  require  customized  transmit  patterns,  implemented  using  only  phase  shifts,  in  order  to 
illuminate  or  suppress  interference  from  specific  regions.  Multiple  custom  patterns  may  be  required  at  every  beam  position,  requiring  either 
precomputation  of  a  large  library  of  phase  weights  or  real-time  computation.  Presented  here  is  a  weighted  pattern-error  fonnulation  that  offers 
efficient  FFT-based  calculation  of  both  the  objective  and  its  gradient  vector  and  is  well  suited  for  solving  with  efficient  quasi-Newton  methods. 
Although  deliberately  simple,  the  objective  has  several  heuristic  parameters  for  flexibility  and  can  be  used  to  simultaneously  customize  the 
mainbeam  shape  and  selectively  suppress  the  sidelobes.  Examples  demonstrate  the  feasibility  of  1000+  element  2D  array  designs,  as  well  as  a 
characterization  of  point-nulls  that  commonly  occur  in  the  mainbeam  of  suboptimal  local  minima. 


15.  SUBJECT  TERMS 

Phase  only  array  pattern 
Array  synthesis 


Transmit  array  pattern 
Array  pattern  optimization 


Phase-only  optimization 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION 

OF  ABSTRACT 

18.  NUMBER 

OF  PAGES 

19a.  NAME  OF  RESPONSIBLE  PERSON 

Dan  P.  Scholnik 

a.  REPORT 

Unclassified 

Unlimited 

b.  ABSTRACT 

Unclassified 

Unlimited 

c.  THIS  PAGE 

Unclassified 

Unlimited 

Unclassified 

Unlimited 

24 

19b.  TELEPHONE  NUMBER  (include  area 
code) 

(202)  404-1943 

Standard  Form  298  (Rev.  8-! 

Prescribed  by  ANSI  Std.  Z39.18 


CONTENTS 


EXECUTIVE  SUMMARY . E-l 

1.  Introduction .  1 

2.  Weighted  Pattern  Error .  2 

2.1  Notation .  2 

2.2  Array  Factor .  2 

2.3  Pattern  Error  Objective .  2 

2.4  Unconstrained  vs.  Constrained  vs.  Feasibility  Formulations .  3 

2.5  Derivation  of  the  Gradient .  3 

2.6  FFT  Computation  of  Objective  and  Gradient .  5 

2.7  Computational  Complexity .  7 

3.  Pattern  Error  Minimization .  7 

3. 1  Local  Solvers .  7 

3.2  Initialization .  7 

4.  Examples .  8 

4. 1  Example  1 :  Circular  Beam .  9 

4.2  Example  2:  Circular  Beam  w/  Sidelobe  Suppression .  10 

4.3  Local  Minima  Characterization .  12 

4.4  Example  3:  Sector  Beam  w/ Nonuniform  Amplitude .  14 

5.  Conclusion .  15 

ACKNOWLEDGMENTS .  16 

REFERENCES .  16 

APPENDIX  A — Implications  of  Symmetry  for  Phase-Only  Array  Factors .  21 


iii 


EXECUTIVE  SUMMARY 


Modern  phased-array  radar  systems  increasingly  require  customized  transmit  patterns,  implemented 
using  only  phase  shifts,  in  order  to  illuminate  or  suppress  interference  from  specific  regions.  Multiple 
custom  patterns  may  be  required  at  every  beam  position,  requiring  either  precomputation  of  a  large  library 
of  phase  weights  or  real-time  computation.  Presented  here  is  a  weighted  pattern-error  formulation  that 
offers  efficient  FFT-based  calculation  of  both  the  objective  and  its  gradient  vector  and  is  well  suited  for 
solving  with  efficient  quasi-Newton  methods.  Although  deliberately  simple,  the  objective  has  several  heuristic 
parameters  for  flexibility  and  can  be  used  to  simultaneously  customize  the  mainbeam  shape  and  selectively 
suppress  the  sidelobes.  Examples  demonstrate  the  feasibility  of  1000+  element  2D  array  designs,  as  well  as  a 
characterization  of  point-nulls  that  commonly  occur  in  the  mainbeam  of  suboptimal  local  minima. 
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A  PARAMETERIZED  PATTERN-ERROR  OBJECTIVE  FOR 
LARGE-SCALE  PHASE-ONLY  ARRAY  PATTERN  DESIGN 


1.  Introduction 

Traditionally,  radar  transmitters  are  run  in  saturation  to  maximize  sensitivity  and  power  efficiency. 
Although  research  into  efficient  amplitude-modulated  transmitters  is  accelerating  [1-3],  in  the  near  term 
high-power  radar  transmitters  will  continue  to  offer  only  phase  control.  This  considerably  complicates  the 
design  of  custom  array  patterns  for  transmit  arrays,  as  all  but  the  most  trivial  phase-only  optimizations 
are  nonconvex  and  thus  may  have  large  numbers  of  local  optima.  In  most  fielded  radar  systems  only  the 
uniformly  illuminated  transmit  beam  is  used,  which  provides  maximum  sensitivity  but  a  narrow  mainlobe 
and  relatively  poor  sidelobes.  In  modem  multifunction/multibeam  arrays,  however,  it  is  often  desirable  to 
trade  sensitivity  for  a  broader  mainbeam,  or  to  selectively  null  directions  or  sectors  of  interest,  or  both.  This 
may  require  patterns  customized  to  each  beam  position,  for  example  to  keep  a  deep  null  around  the  horizon 
as  the  main  beam  is  scanned.  This  in  turn  requires  a  highly  efficient  approach  that  can  be  used  to  precompute 
large  libraries  of  coefficients,  or  to  compute  on  the  fly. 

The  existing  approaches  to  phase-only  array  design  can  be  roughly  divided  into  heuristic,  nonlocal 
optimization,  and  local  optimization  methods.  The  first  two  of  these  are  primarily  of  interest  here  in  that 
they  can  be  used  to  augment  or  initialize  local  minimization.  Many  heuristic  techniques  are  descended 
from  a  well-known  approach  to  nonlinear  FM  waveform  synthesis  [4]  and  only  return  reliably  good  results 
for  impractically  large  arrays  [5].  Examples  of  ID  heuristic  design  approaches  can  be  found  in  [6-11]. 
Reference  [12],  which  proposes  separate  beams  on  interlaced  subarrays,  represents  a  2D  heuristic  method. 
Nonlocal  optimization  methods  (often  referred  to  as  “global”  methods)  such  as  simulated  annealing  and 
genetic  algorithms  have  also  been  used  in  ID  and  2D  [13-18],  although  they  generally  do  not  guarantee  a 
global  solution  in  a  reasonable  time  and  indeed  can  be  impractically  slow  when  applied  directly  to  large 
problems. 

The  drawbacks  of  heuristic  and  nonlocal  methods  suggest  that  local  optimization  should  be  a  part  of 
practical  large-scale  approaches.  Local  optimization  approaches  are  presented  in  [19-24]  for  ID.  The 
reflectarray  and  conformal  literature  is  particularly  rich  with  2D/3D  approaches  [25-29],  many  based  on  the 
alternating -projection  method  of  [30,31].  Many  of  these  approaches  combine  a  core  local  minimization  step 
with  heuristic  or  nonlocal  methods  in  order  to  deal  with  globally  suboptimal  local  minima.  Further  discussion 
on  this  topic  is  deferred  until  3.2.  An  approach  that  seems  to  be  globally  optimal  is  [32],  but  it  is  largely 
restricted  to  custom  sidelobe  nulling  rather  than  mainlobe  shaping.  (An  extension  is  presented  in  [33],  which 
introduces  mainlobe  shaping  but  also  local  minima.) 

Of  the  above  approaches,  local  optimization  combined  with  some  non-local  or  heuristic  methods  to 
avoid  or  discard  “bad”  local  minima  seems  to  be  the  preferred  approach.  To  that  end,  this  work  presents  a 
pattern-error  objective  function  that  mirrors  those  used  in  complex-weight  (non-phase-only)  pattern  design, 
is  flexible,  can  be  efficiently  computed  along  with  its  gradient,  and  is  suitable  for  solution  with  modern  tools 
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that  can  handle  designs  with  thousands  of  elements.  Problems  of  this  class  often  have  local  minima  with 
characteristic  point-nulls  in  the  mainbeam,  which  some  simple  heuristics  arc  used  to  identify  and  avoid. 


2.  Weighted  Pattern  Error 

2.1  Notation 

Throughout  the  sequel,  scalar- valued  variables  (x,<p)  will  be  represented  with  a  medium  weight  italic 
font,  while  vector  and  matrix  values  will  be  represented  in  bold  upright  (x)  or  Greek  (<p)  characters.  The  nth 
element  of  a  vector  x  is  represented  either  explicitly  as  [x]„  or  implicitly  as  xn. 

2.2  Array  Factor 

We  assume  in  the  sequel  that  we  have  identical  isotropic  elements  in  a  planar  array,  so  that  the  array 
pattern  reduces  to  the  2D  array  factor 


N- 1 

A(u)  =  ^  a(xn)eJ^uTxn,  (1) 

n=  0 

where  a(x„)  is  the  complex  array  weight  applied  to  the  element  at  array -plane  location  x„.  Here,  u  =  (  "  j 
represents  two-dimensional  spatial  frequency  normalized  so  that  the  unit  circle  ||u||  <  1  corresponds  to 
propagating  radiation  (so-called  “visible”  space).  The  array  factor  has  the  form  of  a  discrete-space  Fourier 
transform  but  with  an  inverse-transform  exponent  to  describe  plane-wave  propagation. 

The  identical-element  assumption  will  not  always  apply,  and  it  is  straightforward  to  incorporate  unique 
element  patterns  into  the  array  pattern.  When  it  does  apply,  however,  the  benefits  are  substantial:  the  array 
factor  can  be  computed  using  FFT  techniques.  This  enables  efficient  computation  of  key  quantities  used  in 
the  optimization  of  the  array  pattern. 

2.3  Pattern  Error  Objective 

A  very  common  approach  to  the  design  of  FIR  filters,  phased  array  weights,  and  related  structures  is  to 
minimize  some  weighted  error  between  the  response  being  designed  and  some  ideal  desired  response.  When 
coefficient  amplitudes  are  not  fixed,  this  often  results  in  a  convex  optimization  [34, 35].  A  similar  approach 
can  be  used  for  phase-only  responses,  but  will  not  be  convex.  Consider  first  the  weighted  Lp  pattern  error 
metric 


/(a)  =  (  f  W™  (u) 1 1 A(u) |«  -  Dq  (u) f  du)  \p>  1  (2) 

W  '  q  =  {  1,2} 

where  a  is  the  vector  of  array  weights  with  [ a ] „  =  a(x„),  IT  is  a  nonnegative  weighting  function,  and  D  is 
the  desired  magnitude  pattern.  The  support  of  W  determines  the  region  of  integration,  and  its  value  defines 
how  much  each  direction  adds  to  the  overall  error.  It  can  be  zero  for  don't-care  regions.  The  parameters  q 
and  p  permit  some  additional  control  over  the  nature  of  the  resulting  approximation,  with  q  choosing  whether 
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to  approximate  either  the  magnitude  or  power  pattern1  and  real-valued  p  controlling  the  relative  importance 
of  large  and  small  errors.  In  particular,  larger  p  emphasizes  peak  approximation  errors  (thus  providing  more 
of  an  equiripple  error)  while  smaller  p  emphasizes  the  average  errors.  (The  least-squares  special  case  of 
p  =  2  is  already  widely  used  for  phase-only  synthesis,  for  example  [26, 36, 39].)  One  key  difference  from 
traditional  convex  weighted  error  metrics  is  that  only  the  magnitude  of  A  is  considered;  the  phase  response 
will  be  determined  by  the  optimization.  In  complex-weight  design  the  absolute  value  is  usually  avoided  by 
imposing  conjugate  symmetry  on  the  weights,  but  conjugate  symmetry  is  precluded  in  all  but  the  most  trivial 
many  phase-only  designs.  (Details  are  in  Appendix  A.) 

Although  (2)  is  written  in  terms  of  the  array  factor  only,  it  is  straightforward  to  incorporate  a  common 
nonisotropic  element  pattern  Aei  by  defining  VT(u)  =  W'(u)|Aei(u)|  and  D(u)  =  D'lujIAefiujr1.  Thus  the 
simplifying  isotropic  element  assumption  does  not  lead  to  any  significant  loss  of  generality. 

2.4  Unconstrained  vs.  Constrained  vs.  Feasibility  Formulations 

The  focus  here  is  on  unconstrained  minimization  of  (2)  for  absolute  efficiency.  This  is  important  for 
phased-array  applications  as  many  thousands  of  different  beams  might  need  to  be  precomputed,  or  custom 
beams  may  be  required  to  be  computed  on  the  fly.  Another  advantage  of  unconstrained  minimization  is  that  it 
(by  definition)  guarantees  a  feasible  solution.  A  third  advantage  is  that  it  guarantees  (locally)  Pareto-optimal 
solutions;  that  is,  no  part  of  the  pattern  can  be  improved  without  degrading  some  other  part.  This  ensures  that 
all  available  degrees  of  freedom  are  used.  The  primary  drawback  of  using  (2)  for  unconstrained  formulation 
is  that  it  relies  on  the  heuristic  choice  of  a  weighting  function  W  in  order  to  control  the  relative  errors  in 
different  parts  of  of  the  pattern.  Thus  without  some  trial-and-error  or  additional  information  one  cannot 
guarantee  specific  error  bounds. 

A  common  alternative  is  the  family  of  projection-based  feasibility  methods  based  on  [31]  which  have 
no  objective  function  as  such  but  seek  to  directly  meet  a  set  of  constraints.  These  methods  have  their 
own  drawbacks.  First,  the  constraints  may  be  infeasible.  Alternately,  the  constraints  may  be  too  loose,  in 
which  case  the  first  feasible  iterate  found  is  returned  even  though  additional  improvement  is  possible.  In 
general  feasibility  methods  do  not  produce  Pareto-optimal  solutions.  These  methods  are  often  employed  for 
reflectarrays,  where  one  or  a  small  number  of  fixed  beams  are  to  be  designed  and  computational  speed  is  less 
of  a  factor. 

Finally,  there  exist  several  approaches  to  extend  the  unconstrained  optimization  presented  here  to  full 
constrained  optimization  [40].  This  can  be  done  by  adapting  the  weighting  function  W  to  achieve  the  desired 
levels,  which  is  effectively  a  penalty  method.  The  augmented  Lagrangian  method  is  a  related  alternative  that 
can  improve  conditioning  by  estimating  the  Lagrange  multipliers  (dual  variables).  The  drawback  of  these 
methods  is  primarily  increased  computational  complexity. 

2.5  Derivation  of  the  Gradient 

We  now  proceed  to  derive  the  gradient  of  the  objective  /,  which  we  need  for  most  optimization  solvers. 
The  objective  is  a  real- valued  function  of  the  complex  vector  a,  which  we  decompose  as  magnitude/phase 

‘It  has  been  claimed  that  q  =  2  reduces  the  problem  of  undesirable  local  minima,  based  on  results  from  related  analysis 
problems  [36-38].  However,  no  significant  difference  in  local-minima  performance  between  q  -  1  and  q  =  2  was  observed  here. 
This  may  be  attributable  to  the  differences  between  the  analysis  and  synthesis  problems  (including  a  lack  of  a  phase-only  constraint 
in  the  former)  or  to  the  nature  of  the  local  minima  most  often  encountered  in  this  work. 
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p  o  eJI *  where  o  is  the  Hadamard  (pointwise)  product  and  exponentiation  by  a  vector  is  taken  pointwise.  The 
gradient  with  respect  to  both  magnitude  and  phase  can  be  written  as  the  2N  x  1  vector 

V/(p,0)  =  f(p,<f>) 

and  straightforward  application  of  the  chain  rule  to  (2)  yields 
V/  =  fl~P  f  W^(u)||A(u)r  -  D«(u)|P_1sgn(|A(u)|«  -  D^(u))||A(u)|«-2  V|A(u)|2du.  (3) 

It  remains  to  compute  the  gradient  of  |A(u)|2: 

V|A(u)|2  =  v(Re{A(u)}2  +  Im{A(u)}2)  (4) 

=  2Re{A(u)}  VRejA(u)}  +  2Im|A(u)}  VIm|A(u)}. 

Since  A  is  a  holomorphic  (analytic)  function,  it  obeys  the  polar  form  of  the  Cauchy-Riemann  equations: 

p  O  VpRe{A(u)}  =  V0Im{A(u)} 
p  O  Vp  Imj  A(u)}  =  -V^  RejA(u)}. 

Thus  the  gradient  is  completely  defined  by  the  real  and  imaginary  parts  of 

VpA(u)  =  eJ<^  o  eJ^~x  u 

where  X  is  the  matrix  whose  columns  contain  the  in-plane  element  locations  jxo, . . .  ,xjv-i ).  Substituting  back 
into  (4)  and  applying  the  identities  Rc{x  y*\  =  Rej.r)  Rejr/}  +  Imjx}  Imj  y )  and  Im {xy*}  =  Imjx}  Reji/}  - 
Rejx}  Imj  y}  we  find 


(y  \  2  Re{A(u)e  ^  o  e  j~*x'  U1  \ 

\V  J  |A(U)I  =  2  [lm{A(u)p  o  e~J*  °  e~J2fxT"} )  ' 


If  we  now  fix  the  weight  magnitudes  p  we  are  left  with  the  IV  x  1  gradient  | ,4  ( u)  | 2 .  Substituting  back  into 
(3),  we  find  the  nth  gradient  component  as 


df_ 

d(f>n 


=  ImjJ"  f1~pWpq(u)\\A(u)\q  -  Dq(u)\P  1  x 

sgn(|A(u)|<?  -  Dq  (u) )  pne~J(pn  | 


(5a) 
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where  in  (5b)  we  combine  all  the  weighting  factors  into  a  single  function  B  (which  has  the  same  support  as 
W)  and  see  that  the  integral  has  the  form  of  a  Fourier  transform.  Although  this  expression  is  a  little  messy,  it 
is  straightforward  to  compute.  In  practice  the  integrals  in  the  objective  and  gradient  will  be  approximated 
by  a  sum  over  lattice  points.  As  shown  in  the  next  section,  choosing  the  approximating  lattice  properly  in 
relation  to  the  element  lattice  allows  the  forward  and  inverse  transforms  to  be  efficiently  computed  using  an 
FFT/IFFT  pair. 

A  couple  of  concluding  notes:  Although  the  restriction  p  =  1  is  most  common,  representing  an  array  with 
fixed  uniform  amplitude  illumination,  phase-only  optimization  can  also  find  application  to  arrays  with  fixed 
but  nonuniform  tapers.  Such  fixed  tapers  might  be  designed  over  the  ensemble  of  desired  beams,  as  in  [41], 
or  imposed  by  a  legacy  system.  Also,  (2)  can  be  somewhat  trivially  extended  to  a  sum  of  weighted  Lp  norms 
with  different  values  of  p.  This  would  allow  designs  with  different  characteristics  in  different  regions:  an 
equiripple  mainbeam  and  a  L2  sidelobe  region,  for  example.  This  would  still  require  only  a  single  FFT/IFFT 
pair  to  compute. 


2.6  FFT  Computation  of  Objective  and  Gradient 


Because  the  objective  and  gradient  will  be  repeatedly  evaluated  inside  of  an  optimization  routine,  we 
want  to  find  an  efficient  way  to  compute  them.  The  majority  of  the  computation  lies  in  the  sum  of  (1)  and  the 
integral  (to  be  approximated  with  a  sum)  of  (5b).  The  identical-element  assumption  ensures  that  these  are 
both  Fourier  transforms,  and  thus  fast  algorithms  are  available.  If  the  element  locations  are  restricted  to  a 
lattice  (which  they  are  in  most  large  linear  and  planar  arrays),  then  the  computation  can  be  made  using  any  of 
a  number  of  multidimensional  FFT/IFFT  algorithms  [42].  For  arbitrary  element  locations  nonuniform  FFT 
algorithms  exist  [43-45]  that  have  the  same  asymptotic  complexity  as  the  conventional  FFT,  although  the 
overhead  and  scale  factor  arc  greater.  These  find  application  for  conformal  [26]  and  otherwise  irregular  [46] 
arrays,  where  a  uniform  lattice  is  not  available. 

We  focus  attention  here  on  planar  arrays  with  a  uniform  element  lattice;  the  restriction  to  ID  or  extension 
to  3D  is  mathematically  straightforward  although  shadowing  effects  will  generally  violate  the  ideal  element 
pattern  assumption  in  3D.  The  approach  outlined  below  maps  samples  on  an  arbitrary  lattice  such  that  an 
ordinary  rectangular  FFT  can  be  used,  since  such  routines  tend  to  be  widely  available  and  are  the  mostly 
highly  optimized  for  computation.  Similar  approaches  can  be  found  in  [42]. 

Let  2x2  matrix  A  with  linearly  independent  columns  be  a  lattice  basis  for  the  wavelength-normalized 
in-plane  array  element  positions,  such  that  jx„ , . . . ,  xjy-i }  c  A  AZ2  where  Z2  is  the  set  of  all  integer  2-vectors. 
This  lets  us  re-enumerate  the  element  positions  in  terms  of  index  n  e  Z2  as  xn  =  /l  \n.  The  columns  of  the 
dual  lattice  basis  A  7  define  two  adjacent  sides  of  a  parallelogram  that  encloses  one  period  of  the  array  factor. 
We  wish  to  sample  the  array  factor  on  a  uniform  grid  within  this  period  of  sufficient  density  to  approximate 
the  integral  of  (2).  We  do  this  by  defining  a  superlattice  of  the  dual  lattice  with  basis  A~7  R  1 ,  where 
R  =  diag( R] ,  IP)  is  a  diagonal  integer  resampling  matrix.  This  then  defines  |R|  =  R \  R2  spatial  frequency 
points  within  the  period  of  the  form  Uk  =  A~7  R  1  k.  with  k  ranging  over  the  set 


7C  =  jk  =  {kuk2f  :  k]  =  0,...,/?i  -  1;  k2  =0 ,...,R2  -  1). 
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Evaluating  (1)  at  these  points  and  substituting  for  x„  in  terms  of  the  lattice  basis  in  the  exponential  yields 

A(uk)  =  ^  a(xnV2rf  R~‘n 

neZ2 

=  ^  a'(xn)^k'  R"n  (6) 

n  etc 

where 

&  (X„)  =  ^  ^  tl(X(n-Rm))  ( 1 ) 

msZ2 


is  the  periodic  “wrapping”  of  the  aiTay  coefficient  function  a.  For  any  practical  R  the  sum  of  (7)  will  have  at 
most  one  nonzero  term,  and  so  the  operation  (sometimes  also  referred  to  as  data  turning )  represents  only 
zero-padding  and  data  rearrangement.  The  expression  in  (6)  is,  to  within  a  scale  factor,  the  standard  separable 
2D  IDFT,  which  can  be  made  more  explicit  by  breaking  it  out  as 


7f(u(fci,fc2)) 


R  i-l 


Rn  —  l 


z  ^  z 


e  -  a  (X(„  j  n2)) . 


n  i=0 


n  2=0 


This,  to  within  a  possible  normalizing  scale  factor  of  1/|R|,  is  the  form  computed  by  most  standard  2D  IFFT 
routines. 


Having  computed  the  array  factor,  we  now  approximate  the  integral  of  (2)  as  a  Riemann  sum  over  points 
in  the  lattice  A_rR_1Z2: 


/(a)  *  iaHrT  W(Uk)|lA(Uk)l*  -  ^k)f 

where  the  leading  factor  is  the  area  of  the  cell  occupied  by  each  point  uk.  Although  in  the  most-common  case 
W  and  D  would  have  support  over  exactly  one  period  of  the  array  factor,  it  is  not  required.  In  cases  of  wide 
element  spacing,  for  example,  the  visible  space  ||u||  <  1  is  not  contained  within  a  single  period  of  the  array 
factor.  In  such  cases  it  may  make  sense  to  define  the  support  of  W  and  D  over  the  full  visible  part  of  the  array 
factor. 


Similarly  approximating  the  integral  of  (5b)  with  a  sum  yields 


df(</>) 

d(pn 


"’’{WlRI  Z 

V  ksZ2  J 


(8) 


where 


#'(Uk)  =  ^  5(U(k-Rm)) 


(9) 
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is  the  periodic  wrapping  of  B.  (In  this  case  the  sum  of  (9)  may  contain  more  than  one  nonzero  term  if  W 
is  defined  beyond  one  period.)  The  parenthesized  expression  in  (8)  is  the  standard  form  for  the  DFT.  The 
gradient  can  therefore  be  found  by  first  computing  A  via  IFFT  per  (6),  then  computing  B,  and  then  (8)  via 
FFT. 

2.7  Computational  Complexity 

A  straightforward  (non-FFT)  discretization  of  the  gradient  computation  of  (5)  is  essentially  a  matrix- 
vector  multiply,  and  would  have  complexity  O(NM),  where  M  is  the  number  of  points  used  to  approximate 
the  integral.  Using  the  FFT  as  outlined  above  has  complexity  (9(|R|  log2 ( |R| ) ) .  For  most  cases  of  interest 
the  FFT  approach  will  be  considerably  faster,  the  exceptions  being  cases  of  small  N  or  M  that  will  still  be 
quite  fast  regardless  of  the  method.  We  can  make  the  comparison  more  directly  by  noting  that  the  number  of 
sample  points  |R|  is  generally  chosen  proportional  to  the  number  of  elements  N  (multipliers  of  1 00M00  arc 
typical),  so  that  the  asymptotic  complexity  can  also  be  written  simply  as  0(N  log2((V)). 


3.  Pattern  Error  Minimization 

3.1  Local  Solvers 

A  large  number  of  solvers  exist  to  find  local  minima  of  an  objective  function  given  its  gradient.  The 
so-called  quasi-Newton  methods  arc  often  preferred,  as  they  offer  much  better  convergence  properties  than 
simple  steepest-descent  algorithms  but  do  not  require  the  overhead  of  Hessian  matrix  computation  or  inversion 
as  do  true  Newton  methods.  Quasi-Newton  methods  instead  update  an  estimate  of  the  inverse-Hessian  or 
its  equivalent  based  only  on  the  objective  and  its  gradient.  The  original  approaches,  such  as  the  eponymous 
BFGS  method  of  Broyden,  Fletcher,  Goldfarb,  and  Shanno  [47-50]  still  require  storage  of,  and  multiplication 
by,  a  full  N  x  N  matrix,  with  a  per-iteration  complexity  of  0(N2).  These  algorithms  have  previously  been 
used  for  phase-only  array  synthesis  [26,36].  More  recently  limited-memory  versions  (L-BFGS)  [51]  have 
been  developed  that  store  only  a  much  smaller  N  x  L  matrix  to  represent  the  Hessian,  with  per-iteration 
complexity  of  0{NL).  The  penalty  vs.  full  BFGS  is  typically  a  small  increase  in  the  number  of  required 
iterations  to  reach  convergence. 

Several  quasi-Newton  solvers  with  Matlab  interfaces  were  tested  on  the  example  designs  in  the  next 
section.  Several  L-BFGS  implementations  are  freely  available  [52-57].  Of  these,  [54,57]  also  provide 
conjugate-gradient  methods,  an  alternative  class  of  algorithms  for  gradient-based  minimization  with  similar 
computational  complexity.  The  fminunc  function  of  the  commercial  Matlab  optimization  toolbox  provides 
a  full-memory  BFGS  implementation,  although  it  was  found  to  be  substantially  slower  than  the  L-BFGS 
solvers.  All  solvers  showed  similar  convergence  performance,  with  the  L-BFGS  implementations  generally 
requiring  the  least  total  time. 

3.2  Initialization 

As  will  be  seen  in  the  examples,  the  pattern-error  objective  has  some  (indeed  many)  local  minima  that  are 
globally  suboptimal.  A  common  feature  of  these  is  characterized  by  example  in  Sec.  4.3.  Since  the  various 
solvers  find  a  local  minima  “close”  to  the  stalling  point,  choosing  the  initial  element  phases  becomes  key. 
Even  for  simple  designs,  experimentation  quickly  shows  that  choosing  completely  random  phases  rarely 
produces  good  solutions. 
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3.2.1  Approaches  Used  in  the  Examples 

Although  solving  the  local  minima  problem  is  not  the  primary  goal  of  this  paper,  several  heuristic  methods 
were  used  to  improve  local  minima  quality.  These  include: 

Symmetry 

Imposing  even  symmetry  on  the  initial  phases  can  often  improve  local  solutions,  even  if  the  ideal  pattern 
is  not  itself  symmetric.  Circular  symmetry  can  also  prove  useful,  reducing  a  2D  initialization  choice  to  ID. 

Phase  Continuity 

Stalling  with  completely  random  phases  implies  starting  with  a  random  array  factor  which  bears  no 
resemblance  to  the  desired  pattern.  Often  a  better  choice  is  to  sample  a  smooth  phase  function,  such  as  a 
low-order  polynomial.  A  second-order  bivariate  polynomial  will  produce  a  linear- FM  pattern  of  the  form 
considered  in  [5],  with  the  polynomial  coefficients  controlling  the  stalling  beam  position  and  beamwidth. 
(See  Section  4. 1  for  an  example.)  Higher-order  polynomials  provide  a  richer  set  of  initial  phases  but  are  less 
smooth. 

Partial-Convergence  Search 

Often  it  is  apparent  that  the  current  local  solution  is  poor  after  relatively  few  iterations,  well  before  full 
convergence.  One  way  to  exploit  this  is  to  try  many  initializations,  optimize  for  a  fixed  number  of  iterations 
or  until  a  trigger  is  activated  (one  such  trigger  is  discussed  in  Section  4.3),  and  then  fully  optimize  the  best 
partial  result.  These  test  optimizations  can  be  further  sped  up  by  tweaking  parameters  in  the  solvers  (for 
example,  the  memory  depth  L  in  a  L-BFGS  algorithm)  as  well  as  reducing  the  number  of  sample  points  |R| 
that  arc  used  to  approximate  the  integrals  in  the  objective  and  gradient. 

3. 2. 2  Other  approaches 

Although  not  considered  further  here,  more-elaborate  approaches  have  been  suggested  to  mitigate  the 
pattern-synthesis  local  minima  problem,  many  coming  from  the  reflectarray  community.  Some  forego  local 
optimization  altogether,  choosing  instead  metaheuristic  approaches  such  as  simulated  annealing,  particle- 
swarm  optimization,  and  genetic  algorithms  [13-17],  The  primary  problem  with  these  approaches  is 
speed;  [17],  for  example,  describes  a  single  example  that  took  44  hours  to  compute  on  a  single  CPU  core. 
An  approach  in  [16]  that  should  be  applicable  to  local  minimization  as  well  is  to  slowly  walk  the  design 
parameters  out  from  a  uniform-illumination  stalling  point,  re-optimizing  at  each  step. 

Local  optimization  approaches  generally  involve  a  series  of  modified  optimizations.  One  such  approach 
is  to  first  reduce  the  problem  order  by  deriving  the  array  weights  from  a  polynomial  basis  with  optimized 
coefficients  [24,58],  followed  by  full  optimization.  Others  include  adaptively  changing  the  pattern-error 
weighting  [39],  stalling  with  a  significant  amplitude  taper  which  is  gradually  reduced  [28],  or  designing  a 
series  of  increasing-size  arrays  [29].  Both  global  and  local  approaches  are  used  in  [26],  with  the  global  search 
applied  to  low-order  polynomial  models  and  the  result  used  to  seed  the  final  full  optimization. 
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4.  Examples 

In  this  section  a  series  of  examples  is  presented,  with  progressive  complexity.  A  common  array  geometry 
is  used  for  all  examples.  The  array  plane  is  inclined  so  that  boresight  points  15°  above  the  horizon.  The 
basis  matrix  for  the  elements  in  the  plane  is  A  =  (  b  vri  b(2  V3)  ^  which  defines  a  hexagonal  lattice  with  an 

interelement  spacing  of  A/  V3.  The  array  elements  are  centered  on  the  1075  points  of  the  lattice  that  fall 
within  a  radius  of  1 0/1  from  the  origin.  The  1  dB  and  3  dB  beamwidths  of  the  array  with  uniform  (identical) 
drive  phases  arc  about  1.7°  and  3.0°,  respectively. 

For  all  examples,  the  resampling  matrix  of  section  2.6  was  set  to  R  =  diag(512,512),  which  results  in  just 
fewer  than  244  spatial-frequency  sample  points  per  array  element.  Under  Matlab  on  a  dual-core  Intel  Core  i7 
laptop,  the  resulting  objective  computation  took  an  average  of  13  ms,  while  simultaneous  computation  of  the 
objective  and  the  gradient  took  25  ms.  FFT  and  IFFT  computation  contributed  just  under  half  of  the  total 
time. 


4.1  Example  1:  Circular  Beam 


For  a  first  example,  suppose  we  wish  to  form  a  broad  flat -top  circular  boresight  beam  of  radius  10°,  which 
is  approximately  ru  =  0.17  in  units  of  normalized  spatial  frequency.  The  desired  response  D  is  thus  constant 
over  the  circle  and  zero  elsewhere.  By  definition  the  total  energy  in  the  array  weights  is  N,  and  applying 
Parseval’s  relation  reveals  that  the  energy  in  one  period  II  =  {A-  7  v  :  v  e  ([0, 1)  x  [0, 1)) }  of  the  array  factor 
is 


|A(u)|2  du  =  lAl^iV. 


(10) 


If  this  energy  were  distributed  evenly  across  the  circular  beam,  the  resulting  array-factor  height  would  be 
Do  =  J |A|^r:  •  This  is  physically  unrealizable,  of  course,  as  it  implies  zero  sidelobes  everywhere.  Because 
we  know  that  some  of  the  total  energy  must  lie  in  the  sidelobes,  we  will  define  the  relaxed  mainlobe  level 
D  |  =  a  Do,  where  0  <  a  <  1  is  the  relaxation  parameter.  (When  a  is  given  in  dB  units  the  operation 
201og10(a)  is  implied.)  This  relaxed  level  will  then  be  used  in  the  optimization  objective. 

In  conventional  complex  array  weight  optimization  we  would  choose  the  error  weight  function  W  to  be 
nonzero  on  both  mainlobe  and  sidelobe  regions,  weighting  each  to  get  the  desired  balance  between  the  errors. 
Here,  however,  the  choice  of  D\  and  the  implicit  energy  constraint  (10)  combine  to  effectively  constrain  the 
sidelobe  region,  so  that  we  can  set  W  to  unity  within  the  mainlobe  and  zero  elsewhere.  We  argue  this  for 
two  different  cases.  If  we  assume  that  the  optimization  successfully  approximates  the  desired  level  D\  in  the 
mainlobe  region,  then  the  energy  in  the  mainlobe  is  approximately  n r2D2  =  a2\ A|_1  N.  As  a  result  of  (10), 
this  means  that  the  sidelobe  energy  is  approximately  (1  -  ar2)|A|  1  /V.  So  in  this  case  we  have  implicitly  set 
the  total  sidelobe  energy  and  allowed  the  optimization  to  allocate  it  so  as  to  minimize  mainlobe  error.  On  the 
other  hand,  if  D\  is  set  too  high  for  the  optimization  to  approximate  in  the  mainlobe,  this  implies  that  the 
design  has  run  up  against  the  fundamental  Fourier  limitations  of  finite  arrays.  In  this  case  the  optimization 
is  essentially  pushing  upward  on  the  mainlobe  at  optimality,  and  thus  implicitly  is  simultaneously  trying  to 
minimize  sidelobe  energy.  We  can  augment  the  implicit  energy  constraint  with  explicit  sidelobe  weighting  in 
W,  but  will  defer  that  to  the  second  example. 
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Because  the  desired  pattern  is  circularly  symmetric,  the  initial  phases  were  chosen  to  be  circularly 
symmetric  as  well.  A  second-order  polynomial  in  the  element  distance  from  array  center  was  used,  with  the 
polynomial  coefficients  chosen  from  a  uniform  distribution.  Ten  such  polynomials  were  drawn  and  used  to 
generate  ten  prospective  sets  of  stalling  phase  values.  Then,  for  each  of  the  four  combinations  of  p  e  {2,40} 
and  a  €  {OdB, -1.5  dB}  with  D  and  W  as  described  above  and  q  =  1,  the  objective  (2)  was  minimized  for  a 
maximum  of  20  iterations  per  initial  phase  set.  The  best  result  of  the  ten  for  each  subdesign  was  then  fully 
optimized.  The  time  to  test  all  ten  initializations  using  [52]  was  about  2  s  per  subdesign,  while  the  time  to 
fully  optimize  the  best  candidate  ranged  from  0. 1  s  to  0.6  s.  Convergence  was  typically  achieve  in  less  than 
50  iterations. 

A  special  case  of  the  second-order  polynomial  used  is  a  pure  quadratic  in  the  element  distance.  In  [5]  this 
was  shown  to  amount  to  2D  linear  FM  modulation,  which  in  the  limit  of  an  infinite  circular  array  approaches 
an  ideal  flat-top  circular  beam.  The  single  polynomial  coefficient  sets  the  beamwidth,  and  for  practical  array 
sizes  also  affects  the  mainbeam  ripple  characteristic.  For  comparison  to  the  Lp -optimized  designs,  this  single 
parameter  was  optimized  to  maximize  the  minimum  gain  within  the  mainlobe. 

The  resulting  array  factors  are  shown  in  Fig.  la.  They  are  scaled  relative  to  the  pattern  of  an  isotropic 
antenna  radiating  the  same  power;  for  unit-magnitude  coefficients  this  is  equal  to  A(u)/  VfV-  Visually  all 
four  Lp  designs  look  very  similar  in  Fig.  la,  which  shows  the  patterns  over  all  visible  spatial  frequencies, 
with  differences  primarily  in  the  mainlobe.  The  1  inear  FM  pattern,  by  comparison,  has  significantly  more 
ripple  and  wider  transition  bands.  The  mainlobe  detail  is  shown  in  Fig.  lb,  and  here  we  can  most  clearly 
see  the  effect  of  the  values  of  p  and  a.  When  a  =  1  (OdB),  we  see  that  the  optimized  mainlobe  falls  well 
below  the  target  Do  level  except  at  the  center.  When  we  relax  the  desired  mainlobe  by  1 .5  dB,  then  we  see 
the  optimized  mainlobe  ripple  about  the  desired  level  as  in  classic  digital  filter  designs.  As  we  might  expect, 
the  p  =  2,  a  =  OdB  case  provides  the  greatest  rms  value  in  the  mainlobe,  only  1.1  dB  below  the  ideal  Do 
value,  and  thus  concentrates  the  most  energy  in  the  mainlobe.  However,  it  has  a  worse  minimum  value  in  the 
mainlobe  than  either  p  =  40  case,  which  are  essentially  equiripple  owing  to  the  large  p.  In  many  broad-beam 
applications  it  is  this  maximin  pattern  value,  representing  worst-case  gain,  that  is  most  important.  (Despite 
being  optimized  specifically  for  this  metric,  the  linear  FM  weighting  is  significantly  worse  than  all  of  the 
Lp  designs  on  this  score.)  Adding  -1.5  dB  relaxation  to  the  p  =  40  case  has  minimal  effect  on  either  rms  or 
minimum  mainlobe  values  and  simply  causes  the  energy  from  the  large  central  ripple  to  be  distributed  to  the 
sidelobes.  Further  relaxing  the  mainbeam  height  does  result  in  an  increasingly  flat  mainlobe  at  the  expense  of 
greater  sidelobe  energy. 

4.2  Example  2:  Circular  Beam  w/  Sidelobe  Suppression 

Building  on  the  previous  example,  we  now  enforce  a  zone  of  increased  sidelobe  suppression  from  2° 
below  to  2°  above  the  horizon.  Such  a  zone  could  be  useful  to  mitigate  multipath,  for  example,  or  to  reduce 
interference  with  other  systems.  The  extra  suppression  is  achieved  by  setting  D  to  zero  in  this  region,  and 
setting  W  to  a  nonzero  constant.  Adjusting  the  relative  value  of  W  on  the  mainlobe  and  sidelobe  regions 
allows  us  to  tradeoff  the  residual  errors  in  the  two  regions.  As  before,  the  rest  of  the  sidelobes  are  allowed  to 
be  determined  by  the  optimization.  The  initialization  is  handled  as  in  the  previous  example,  with  parameters 
p  =  2  and  a  =  -1.5  dB.  Although  the  computational  complexity  of  each  iteration  is  identical  to  the  previous 
example,  more  iterations  (100  was  typical)  are  required  for  convergence,  resulting  in  an  average  total  solution 
time  of  3  s.  In  general,  such  closely  spaced  mainlobe  and  sidelobe  specifications  appear  to  produce  error 
surfaces  that  are  not  locally  well  approximated  as  a  quadratic  form  and  result  in  slower  convergence. 
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(a)  Full  visible  space  array  factors.  The  lines  on  the  top  plot  indicate  constant  azimuth  and  elevation  at  15°  intervals,  with  array 
boresight  at  15°  elevation.  The  bottom  row  shows  azimuth  slices  through  the  center  of  the  mainlobe. 


(b)  Mainlobe  detail.  The  lines  on  the  top  plot  indicate  constant  azimuth  and  elevation  at  5°  intervals.  The  bottom  row  shows  azimuth 
slices  through  the  center  of  the  mainlobe.  The  green  (upper)  line  segment  indicates  the  theoretical  ideal  array-factor  height  Dq,  while 
the  red  line  indicates  the  relaxed  height  D\ . 


Fig.  1 :  Optimized  array  factors  for  the  Example  1  designs. 


The  optimized  array  factor  is  shown  in  Fig  2,  with  azimuth  and  elevation  cuts  shown  above  and  to  the 
right.  The  deep  suppression  around  the  horizon  is  clearly  visible,  with  the  peak  and  rms  horizon  sidelobe 
levels  indicated  on  the  plot.  The  max,  min,  and  rms  values  in  the  mainlobe  are  also  shown  on  the  plot,  and 
reveal  that  the  “cost”  of  the  extra  sidelobe  suppression  is  a  loss  of  0.2  dB  of  mainlobe  energy  and  0.7  dB 
of  worst-case  gain.  Visually,  we  can  see  that  the  sidelobes  at  elevations  above  the  mainlobe  have  increased 
relative  to  the  previous  design. 


Figure  3  is  a  different  view  of  the  same  optimized  array  factor,  to  illustrate  the  relationship  with  the 
functions  D  and  W.  The  parallelogram  bounds  one  period  of  the  array  factor,  over  which  W  and  D  are 
defined.  The  two  sides  meeting  at  the  origin  represent  the  columns  of  the  dual  lattice  basis  T  =  (  ^  ° - 
Superimposed  on  the  plot  are  (tiny)  white  dots  at  every  nonzero  sample  of  W. 
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Azimuth  cut  at  v  =  0.00 
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Fig.  2:  Optimized  array  factor  for  the  second  example.  The  plots  above  and  to  the  right  are  azimuth  and 
elevation  slices  through  the  beam  center.  The  red  line  segments  indicate  the  desired  array  factor  level  D\ . 


4.3  Local  Minima  Characterization 

The  complicated  error  surface  for  the  objective  (2)  means  that  numerous  (globally)  suboptimal  local 
minima  exist  for  most  designs.  Local  solvers  will  find  the  “closest”  (in  some  sense)  minimum  to  the  point 
used  to  initialize  the  optimization.  Thus  the  results  can  be  very  sensitive  to  a  change  in  the  initial  phases.  In 
the  examples  presented  so  far,  the  combination  of  a  smooth  and  circularly  symmetric  phase  initialization 
function  and  running  10  partial  trials  was  almost  always  sufficient  to  produce  a  good  local  (and  possibly 
global)  solution.  With  more-complicated  desired  geometries  it  can  be  hard  to  avoid  “bad”  local  minima. 

Here  we  define  “bad”  local  minima  specifically  as  those  with  point  nulls  within  the  mainbeam.  Reference 
[29]  mentions,  somewhat  in  passing,  a  local  minima  characteristic  consisting  of  “holes”  or  point  nulls  in  the 
mainbeam  of  a  pattern  designed  using  a  least-squares  approach.  The  same  nulls  were  seen  both  with  and 
without  the  phase-only  constraint,  and  thus  are  due  to  optimizing  only  the  magnitude  (rather  than  complex) 
pattern.  Point-null  artifacts  have  also  been  described  as  the  Fourier-domain  counterpart  to  the  image-domain 
“stripes”  observed  and  characterized  in  phase-retrieval  problems  [59]. 

Consider  Fig.  4a,  which  results  from  choosing  the  worst  of  the  10  initial  phases  rather  than  the  best  in 
Example  2.  Eight  point  nulls  are  clearly  visible.  (The  indicated  minimum  mainbeam  value  is  an  artifact  of  the 
grid  used  for  plotting  and  computation;  each  point  is  a  true  null  to  working  precision.)  The  mainlobe  details 
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Fig.  3:  Illustration  of  the  geometry  of  the  D  and  W  functions.  The  parallelogram  bounds  one  period  of 
the  array  factor.  Overlaid  onto  the  period  are  (tiny)  white  dots  for  each  nonzero  sample  of  W,  defining  the 
mainlobe  and  horizon  regions. 


in  magnitude  and  phase  are  shown  in  Fig.  4b.  In  the  immediate  orbit  of  each  null,  the  phase  of  the  array 
factor  rotates  a  full  360  degrees  [59].  Along  any  line  drawn  through  a  null  location  the  phase  is  discontinuous 
(at  the  null),  and  since  the  complex  array  factor  (as  a  finite  discrete-space  Fourier  transform)  is  continuous, 
the  magnitude  must  go  to  zero  at  that  point. 


As  noted  in  [29, 59],  the  nulls  are  narrow  enough  that  they  do  not  greatly  affect  a  squared-error  objective. 
One  might  expect  here  that  increasing  p  would  directly  mitigate  the  point  nulls,  since  they  will  then 
increasingly  dominate  the  objective.  Experimentation  has  shown  otherwise,  with  point  nulls  appearing  for 
any  p.  Similar  nulls  arise  when  using  alternating -projection  approaches  [36,59].  Thus  point  nulls  appear  to 
be  a  common  feature  of  magnitude/power  pattern  synthesis  problems,  and  not  specific  to  a  given  formulation. 


There  are  (at  least)  three  exploitable  properties  of  such  nulls:  they  tend  to  form  early  in  the  solver  iteration, 
they  are  stable  once  formed,  and  they  are  easy  to  detect  while  computing  (2).  This  can  be  used  to  trigger  an 
early  solver  termination,  in  order  to  move  on  to  another  starting  point.  For  example,  the  minimum  value  of 
the  ratio  |A(u)|/D(u)  over  the  mainbeam  region  after  20-60  iterations  has  proved  to  be  a  reliable  indicator 
of  point  nulls  [60]. 
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(a)  All  visible  space. 


(b)  Details  of  the  magnitude  and  phase  response  in  the  mainlobe.  The  nulls 
in  magnitude  correspond  to  discontinuities  in  phase. 


Fig.  4:  Array  factor  with  point  nulls  for  a  representative  “bad”  local  minimum  of  the  Example  2  design. 


4.4  Example  3:  Sector  Beam  w/  Nonuniform  Amplitude 

Finally,  we  consider  an  example  with  a  noncircular  mainbeam  which  is  not  centered  at  boresight.  The 
mainbeam  region  is  the  sector  between  0°  and  30°  in  azimuth  and  between  5°  and  60°  in  elevation.  (A  similar 
design  can  be  found  in  [27].)  Rather  than  a  flat  desired  pattern,  we  seek  a  pattern  that  provides  uniform 
radar  detection  sensitivity  at  the  lesser  of  100  km  and  the  range  at  40  km  altitude.  (The  actual  units  are 
irrelevant,  only  the  ratio  matters.)  The  result  is  a  beam  that  is  flat  for  low  elevation  angles  and  falls  off  (as 
cosecant-squared  in  magnitude)  for  higher  elevation  angles  as  the  altitude  limit  kicks  in.  As  before  we  scale 
down  the  desired  response  from  the  ideal  to  account  for  sidelobe  energy,  in  this  case  by  1  dB.  We  retain 
the  horizon  sidelobe  suppression  region  from  Example  2.  Initialization  is  done  similarly  to  the  previous 
examples,  even  though  now  the  desired  mainbeam  is  far  from  circularly  symmetric.  Two  changes  to  increase 
initialization  diversity:  the  polynomial  order  of  the  phase  initializing  functions  is  increased  to  four,  and  50 
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starting  points  are  tested  at  a  time.  Even  with  these  extra  initial  values,  multiple  runs  are  sometimes  required 
to  find  a  single  null-free  solution.  Here  initialization  took  about  9  s,  and  final  solve  took  over  200  iterations 
and  about  1.8  s. 

The  array  factor  for  one  solution  is  shown  in  Fig.  5.  The  esc2  profile  is  clearly  visible  in  the  elevation 
slice,  and  we  achieve  a  deeper  horizon  depth  than  in  Example  2  (at  the  expense  of  higher  sidelobes)  by 
increasing  the  weight  on  that  region  of  W. 


Azimuth  cut  at  v  =  0.00 


dBi 


Fig.  5:  Optimized  array  factor  for  the  third  example.  The  plots  above  and  to  the  right  are  azimuth  and 
elevation  slices  through  the  beam  center. 


5.  Conclusion 

Phase-only  array  pattern  optimization  is  both  less  common  and  harder  to  solve  than  conventional  receive- 
pattem  synthesis.  Modem  tools  have  made  general-purpose  optimization  far  more  accessible  and  efficient, 
and  are  exploited  here  with  an  efficient  yet  reasonably  flexible  phase-only  objective  function  that  mirrors 
common  aproaches  from  the  non-phase-only  world.  The  examples  demonstrate  that  it  is  now  quite  feasible  to 
quickly  and  directly  optimize  thousands  of  element  phases  to  approximate  arbitrary  responses.  The  nonconvex 
pattern-error  objective  has  many  local  minima  with  catastrophic  features  such  as  nulls  in  the  midst  of  the 
mainlobe.  For  simpler  designs  fairly  simple  heuristics  were  demonstrated  to  avoid  such  bad  local  minima. 
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These  were  found  to  be  less  effective  for  more  complicated  designs,  and  so  more  work  is  needed  to  produce  a 
fully  unsupervised  algorithm. 
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Appendix  A 


IMPLICATIONS  OF  SYMMETRY  FOR  PHASE-ONLY  ARRAY  FACTORS 

Symmetries  arc  often  imposed  on  filter  and  array  coefficients,  primarily  for  the  purposes  of  reducing  the 
number  of  optimization  variables  without  significantly  impacting  the  achievable  magnitude  response.  Several 
common  symmetries  and  their  corresponding  restrictions  on  the  array  factors  are  summarized  below: 


a(x)  =  ±a(-x) 

4 — 4  A(u)  =  ±A(-u) 

Even/Odd  symmetry 

Even/Odd  symmetry 

a(x)  =  ±a*(-x) 

4 — >  A(u)  -  ±A*(u) 

Conjugate  symmetry/antisymmetry 

Real/Imaginary 

a(x)  =  ±a*(x) 

4 — 4  A(u)  =  ±A*(-u) 

Real/Imaginary 

Conjugate  symmetry/antisymmetry 

Other  symmetries  are  possible  in  2D,  but  these  will  suffice  for  the  present  discussion.  Of  these,  some  have 
obvious  drawbacks.  A  restriction  to  purely  real  or  imaginary  weights  would  be  of  limited  use  in  the  phase-only 
case,  since  it  would  restrict  the  weights  to  +1  or  ± /,  respectively.  Odd  symmetry  implies  A(0)  =  0,  and  thus 
precludes  any  mainbeam  that  included  boresight. 

Conjugate  (linear-phase)  symmetry  initially  appears  very  attractive  because  it  reduces  the  number 
of  optimization  variables  by  half  with  no  immediately  obvious  restriction  on  the  magnitude  response. 
Indeed,  conjugate  symmetry  is  widely  used  when  array/filter  weights  are  not  magnitude  constrained.  As 
shown  in  the  following,  however,  conjugate  symmetry  combined  with  the  phase  only  constraint  imposes 
a  rather  strict  implicit  constraint  on  the  array  factor  mainlobe.  We  write  the  phase  only  constraint  as 
|a(x)|2  =  a(x)a*(x)  =  I  ^(x),  where  1  #  is  the  indicator  function  on  the  array  support  X.  To  simplify  the 
argument,  we  will  assume  that  X  lies  on  a  lattice  so  that  A  is  periodic.  A  Fourier  convolution  property  then 
allows  us  to  write 


(A  *  At)(u)  =  Ai(u)  (Al) 

where  *  indicates  (periodic)  convolution,  A' (u)  =  A*(-u)  is  the  paraconjugate  operation,  and  A]  is  the  array 
factor  resulting  from  all  unity  weights.  Now  A  is  real  and  thus  entirely  positive  or  negative  valued  on  any 
null-free  mainlobe  region.  Thus  for  a  conventional  array  factor  containing  a  single  mainlobe  surrounded  by 
significantly  lower  sidelobes,  the  left  side  of  (Al)  will  be  a  boresight  beam  with  a  single  mainlobe  of  roughly 
twice  the  width  of  the  mainlobe  of  A.  However,  we  know  that  Ai  has  the  narrowest  mainlobe  of  all  array 
factors  on  the  support  A,  in  the  sense  that  it  maximizes  |A(0)|  for  a  fixed  total  weight  energy  ||a||2.  Thus, 
evidently  imposing  conjugate-symmetric  symmetry  on  phase-only  weights  precludes  synthesizing  a  mainlobe 
significantly  broader  than  that  of  A\.  A  similar  argument  holds  for  conjugate  antisymmetry. 

This  would  appear  to  leave  even  symmetry  as  the  most  viable  symmetry  out  of  those  listed  for  phase-only 
array  factors,  although  its  utility  is  limited  by  the  corresponding  restriction  on  the  array  factor. 
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